Estimating Access to Rabies Post-Exposure Prophylaxis in Madagascar
INTRODUCTION
Access to healthcare is a significant driver of disease burden globally, and heterogeneities in access to care due to socioeconomic and geographic factors likely shape where and who are the most impacted by a disease, particularly at sub-national levels. However, these factors are rarely incorporated into estimates of burden due to limited data.
Deaths due to canine mediated rabies, estimated to cause approximately 60,000 human deaths anually, can be prevented through prompt administration of post-exposure prophylaxis. However, access to this intervention is highly limited in areas where the disease is endemic, often due to inaccessibility to health centres which provide PEP (cite GAVI paper or Nandini’s paper of geographic availability of PEP). Data on true rabies exposures in humans and incidence in animals is also lacking in most of these countries, with the most commonly available data being numbers of bite victims reporting to health facilities. The majority of rabies burden studies use these data to estimate burden of rabies from probability decision tree frameworks, often with the key assumption that overall reported bite incidence (i.e. both bites due to non-rabid and rabid animals) are proportional to rabies incidence (i.e. the morebites reported in a location, the higher the incidence of rabies exposures there) and that reporting is uniform across space. While at a national level, these estimates may be accurate, at the sub-national level, this framework will likely underestimate rabies deaths in places with low reporting and overestimates rabies deaths in places with high reporting of bites. The most recent estimation of burden and the impact of PEP used another approach, using transmission dynamic models as the backbone to predict incidence based on the level of vaccination coverage and size of the dog population at the national level. Using transmission dynamic models to estimate incidence could improve upon previous studies which may underestimate rabies burden in areas with low reporting.
In Madagascar, Institut Pasteur de Madagascar (IPM) provides PEP at no-cost to patients at 31 clinics across the country. PEP is not available at any other public clinics or through the private sector. Due to the spatially restricted nature of PEP, geographic access is likely to be a major driver of disease burden within the country. To get spatial estimates of disease burden in this context, we flip the standard decision tree and make the assumption that reported bite incidence reflects access and reporting to PEP rather than differences in rabies incidence, using travel times to clinics as a predictor of reported bites. Then using a range of rabies incidence given endemic transmission with no mass dog vaccination (GAVI paper), we generate sub-national estimates of rabies burden in an adapted decision tree framework. Finally, using this same model pipeline, we explore the impacts of geographically expanding access to PEP in Madagascar on reducing human rabies deaths.
METHODS
Bite patient data
We used a database of bite patient forms submitted to IPM from 24 ARMC across the country between 2014 - 2017. These were individual patient data forms that were submitted as frequently as monthly to annually by clinics. Two clinics, the IPM ARMC and the Fianarantsoa ARMC used computer databases from which the data during this period were extracted and merged to the larger database. These data include the administrative district of the bite patient’s address and the date of reporting. As previous analyses showed likely substantial undersubmission of forms from clinics (cite baseline data paper), we calculated the proportion of days of the year for which forms were submitted, excluding any periods of time for which there were no forms submitted for 15 consecutive days. Clinic level reporting was calculated the proportion of days for which were included based on the 15 day threshold per year. As low risk contacts with probable or confirmed rabies cases have been shown to make up approximately 20% of patients reporting to ARMC (Rajeev et al. 2018), we attempted to exclude these by looking at the distribution of patients reporting per day to a given clinic. Generally, contacts present as clustered cases, so we excluded patients reporting on any dates that had greater than 3 standard deviations above the mean number of patients reporting per day (Figure 3A). We validated this method using 27 months of bite patient data from the Moramanga ARMC for which we had details on the type of exposure.
GIS Data
We use the global friction surface for 2015 generated by the Malaria Atlas Project ( https://map.ox.ac.uk/research-project/accessibility_to_cities/, Weiss et al. 2015,) and GPS points of clinics to estimate the travel time to the nearest ARMC for the country of Madagascar at a 1 x 1 km scale. We then calculated a weighted average of travel times by human population to the commune level, using administrative shapefiles available trhough the UN Office for the Coordination of Humanitarian Affairs. For each clinic, we defined the catchment area as all districts for which the clinic was the closest ARMC. Human population estimates were taken from the 2015 UN adjusted population projections from World Pop (www.worldpop.org, Linaird et al. 2012) and also aggreagated to the commune level
Model of reported bites as a function of travel time
While the national bite patient data is available at the district level, travel times vary significantly within a district (see Figure SX, need to think about how to show this). In order to translate the impacts of differences in access at sub-district scales to the magnitude of reported bites at the district scale, we modeled bites at the district level as the sum of incidence at the commune level. Incidence at the commune level is then a function of travel times to the closest ARMC. Specifically, we modeled bites as follows: \[ \mu_{d} = \sum \limits_{j=1}^jexp(\beta_{t}T_j + \beta_0)\times pop_j \] where \(\mu_{d}\) is the mean number of bites in district, which is the sum of bites at the commune level given commune level travel times, \(T_j\). We then estimate the likelihood of observing the bites at the district level where bites are a poisson distribution around the mean \(\mu_{d}\), given \(\beta_t\), the effect of travel times of reported bites and \(\beta_0\), the model intercept.
As we had data available on reported bites at the commune level from the Moramanga ARMC (from Rajeev et al. 2018), we modeled observed commune bites in the same framework, but un-aggregated, where: \[ \mu_{j} = exp(\beta_tT_j + \beta_0)\times pop_j \] where \(\mu_{j}\) = the mean number of bites in commune \(j\) and the observed bites at the commune level follow a poisson distribution around the mean \(\mu_d\). We only included communes which were designated to be within the catchment for the clinic. We compared our estimates of \(\beta_t\) (i.e. the impact of travel time on incidence) and \(\beta_0\) (the intercept) for our district data at the national level and the commune level data from the Moramanga ARMC.
Model sensitivity/validation
- Out-of-fit check for Tana + Fianarantsoa?
- Compare to null models
- Look at sensitivity to inclusion criteria of days and exclusion of contacts
Estimation of burden and reporting
We used our model to predict average annual bite incidence for all 114 districts in Madagascar, and estimated average reporting of rabid exposures and deaths due to rabies given this and assumptions about rabies exposures.
We calculated the expected reporting of rabid exposures (\(\overline\rho\)) given bite incidence as predicted by our model(\(\mu\)) $: \[
\overline\rho = \frac{\mu \times p_{rabid}}{R}
\] or the fraction of incidence due to rabid exposures (\(\mu \times p_{rabid}\)) divided by the total rabies exposure incidence (\(r\)) for a range of estimated rabies incidence and \(p_{rabid}\). We look at the range of \(p_{rabid}\) reported in Rajeev et al. 2018 for data from the Moramanga District (0.2 - 0.6). where the proportion of reported bites that are rabies exposures (\(p_{rabid}\)) are defined as: \[
p_{rabid}=
\begin{cases}
x, & \text{if}\ \frac{R \times \rho_{max}}{B} > x \\
\frac{R \times \rho_{max}}{B}, & \text{otherwise}
\end{cases}
\]
such that rabid reported bites (i.e. \(p_{rabid} \times B_i\)) cannot exceed the expected number of human exposures given maximum reporting (i.e. \(R_i \times \rho_{max}\)). \(\rho_{max}\) taken from the Moramanga ARMC data for Moramanga Ville, the commune with the ARMC (i.e. the area with the minimum travel time). .
The rabies incidence in dogs in the absence of any vaccination, \(r\), is multiplied by the estimated dog population in the commune (\(D_i\)) and the exposure rate per rabid dogs (\(p_{exp}\) = 0.39 persons exposed per rabid dog)(Hampson et al. 2018) to generate \(R\). We estimate the dog population by using a human:dog ratio of 5 to generate our maximum expected incidence and an HDR of 25 for our minimum expected incidence. As there is little data on dog populations in Madagascar, this range of HDRs encompasses a wide range observed across Africa (cite!).
To estimate the average number of deaths for each commune, we extended the above framework into a stochastic framework as follows: \[ deaths_i = (R_i - p_{rabid}B_i) \times p_{death} \] where \(R_i\) is drawn from a uniform distribution between the minimum and maximum expected number of human exposures for each location and \(B_i\), the number of reported bites, is drawn from a poisson distribution with the mean predicted number of bites from the travel time model. We assume that all rabies exposures reported to an ARMC receive and complete PEP, and PEP is completely effective at preventing death due to rabies.
Need to include bits about RIG in here somewhere
Estimating the impact of expanding PEP Access
We use this framework to compare three scenarios of PEP provisioning in Madagascar:
- The baseline with the current clinic locations (n = 31)
- Expansion of ARMC to one clinic per district (n = 114)
- Expansion of ARMC to all CSB IIs (something to describe CSB IIs here), that is primary hospitals (n = X)
We use data on the location of CSBs provided by IPM to regenerate travel times to the nearest ARMC given expansion as per scenario 2 and 3. We then predict the expected bite incidence from the model given these new travel times and compare the relative decreases in burden for the three scenarios.
Sensitivity analyses for burden estimates
- R (increasing or decreasing with pop) (min + max)
- P-rabid min and max
- Probability of death min and max
- Look at change in raw # of deaths and also change in burden by scenarios
RESULTS
Estimating average bite incidence at the district level
Figure 1C
Figure 1B
Figure 1C
For most districts, the majority of bites were reported to the closest clinic as estimated by our travel time metric (Fig 1B), although some patients did report to the ARMC that were not closest to them in terms of travel times (Figure 1A). In addition, at the clinic level, the majority of bites reported to a given ARMC came from within the catchment (Figure 1C).
Figure 2A
Figure 2B
We estimated clinic level reporting as the proportion of days on which forms were submitted, excluding any periods for which there were no forms submitted for 15 consecutive days (Figure 2A). This estimate did vary based on our assumption of the threshold number of consecutive days (Figure 2B), so we looked at sensitivity of parameter estimates to changing this threshold (Supplementary Materials?). In all cases, we further excluded any years for which there was less than 25% reporting.
Figure 3A
Figure 3B
We attempted to exclude these by looking at the distribution of patients reporting per day to a given clini, excluding patients reporting on any dates that had greater than 3 standard deviations above the mean number of patients reporting per day (Figure 3A). We validated this using the Moramanga ARMC data for which we had details on the type of exposure, and found that setting the threshold to 3 standard deviations resulted in approximately 50% of known contacts excluded, with only 2% of non-contacts excluded (Figure 3B). For the national data for which a subset of patient forms were noted to be contacts, we found that our exclusion criteria of 3 SDs identified 68.28 of known contacts. We further excluded known contacts which were not identified based on the daily distribution of patients, resulting in the exclusion of approximately 7.18% of patient records from the national data.
Figure 4
Our final dataset consisted of estimates of average bite incidence for 75 districts from 21 catchments (Figure 4).
Models of bites as a function of travel times
Figure 5A
| model | par | values | upper_CI | lower_CI |
|---|---|---|---|---|
| Mada | B_ttimes | -0.0099980 | -0.0098088 | -0.0101782 |
| Mora | B_ttimes | -0.0125278 | -0.0114495 | -0.0137090 |
| Mada | intercept | -4.9853776 | -4.9683780 | -5.0038438 |
| Mora | intercept | -5.2735226 | -5.1597086 | -5.3831747 |
Figure 5B
We estimated similar parameter values from our commune-level data from the Moramanga ARMC and the district level data from 19 clinics across the country (Table 1, Figure 5A), with reported bite incidence decreasing with travel times to the ARMC. Both models produced reasonable fits to the data (Figure 5B,…), however, there was some variation in bite incidence which both models were not able to capture.
Model sensitivity analyses
- compare to null models: flat incidence, incidence solely predicted by pop;
- sensitivity of param estimates to reporting cut-offs (i.e. days excluded)
- sensitivity of param estimates to exclusion of contacts
Estimation of burden and reporting
Figure 6
Generally, estimated reporting of rabies exposures decayed with travel times given model predicted bite incidence and a range of rabies incidence and \(p_{rabid}\) (Figure 6). Given our model assumptions, reporting was estimated at the maximum of 0.98 for travel times under 1 hour given the maximum estimated rabies exposure incidence and the minimum estimate of \(p_{rabid}\) (the lower range of reporting probabilities), and travel times under 5.5 hours given the minimum estimated rabies exposure incidence and the maximum estimate of \(p_{rabid}\) (the upper range of reporting probabilities).
Figure 7
When we estimate burden of deaths stochastically within this range of incidence and given our high and low estimates of proportion of reported bites that are rabid, we see that burden of deaths also decreases with travel times (Figure 7, results presented aggregated at the district level). Overall, we estimate average annual deaths between 405 - 733. Also estimate deaths averted here!.
Burden sensitivity analyses
- Compare to assuming flat reporting + bite model? Burden @ district level with bites @ district level?
- \(\rho_{max}\), change by urban vs. rural commune?
- If rabies incidence, \(R\), scales with population size within our estimated range
- Impacts on # of deaths vs. change in scenarios (i.e. focusing on relative decreases in incidence of deaths based on expanded access rather than raw #s)
Discussion
Key findings
Strengths and Limitations
- Not accounting for clinic functioning
- Not explicit data on reporting or rabies incidence!
- Other factors that drive reporting
- While we use a signficiant amount of assumptions, these are the same set of assumptions that underlie current rabies estimates, we’re just improving them a bit… and ours are a likely a little bit more in line with reality @ a sub-national level?